Para esta sección seguiremos el libro Doing Bayesian Data Analysis de John K. Kruschke.
Hasta ahora hemos estudiado métodos estadísticos frecuentistas (o clásicos). El punto de vista frecuentista se basa en los siguientes puntos (Wasserman, 2005):
La probabilidad se refiere a un límite de frecuencias relativas, las probabilidades son propiedades objetivas en el mundo real.
En un modelo, los parámetros son constantes fijas (desconocidas). Como consecuencia de que los parámetros son constantes no se pueden realizar afirmaciones probabilísticas en relación a éstos.
Los procedimientos estadísticos deben diseñarse con el objetivo de tener propiedades frecuentistas bien definidas. Por ejemplo, un intervalo de confianza del 95% debe contener el verdadero valor del parámetro con frecuencia límite del 95%.
Por su parte, el paradigma Bayesiano se basa en los siguientes postulados:
La probabilidad describe grados de creencia, no frecuencias limite. Como tal, uno puede hacer afirmaciones probabilísticas acerca de muchas cosas y no solo acerca de datos sujetos a variabilidad aleatoria. Por ejemplo, puedo decir: “La probabilidad de que Einstein tomara una copa de te el 1ro. de agosto de 1948 es del 35%”; esta proposición no hace referencia a ninguna frecuencia relativa sino que refleja la certeza que yo tengo de que sea verdadera.
Podemos hacer afirmaciones probabilísticas de parámetros.
Podemos hacer inferencia estadística de un parámetro \(\theta\) por medio de distribuciones de probabilidad. Las inferencias como estimaciones puntuales y estimaciones por intervalos se pueden extraer de dicha distribución.
Kruschke describe los puntos de arriba como dos ideas fundamentales del análisis bayesiano:
La inferencia bayesiana es la reubicación de creencias a lo largo de posibilidades.
Las posibilidades son valores de los parámetros en modelos descriptivos.
¿Qué tanta certeza tienes de que una moneda acuñada por la casa de moneda mexicana es justa? Si, en cambio, consideramos una moneda antigua y asimétrica, ¿creemos que es justa? En estos escenarios no estamos considerando la verdadera probabilidad, inherente a la moneda. Lo que queremos medir es el grado en que creemos que cada probabilidad puede ocurrir, en lugar de una serie de experimentos llevados a cabo con las monedas.
Para especificar nuestras creencias debemos medir que tan verosímil pensamos que es cada posible resultado. Describir con presición nuestras creencias puede ser una tarea difícil, por lo que exploraremos como calibrar las creencias subjetivas.
Considere una pregunta sencilla que puede afectar a un viajero: ¿Qué tanto crees que habrá una tormenta que ocasionará el cierre de la autopista México-Acapulco en el puente del 20 de noviembre? Como respuesta debes dar un número entre 0 y 1 que refleje tus creencias. Una manera de seleccionar dicho número es calibrar las creencias en relación a otros eventos cuyas probabilidades sean claras.
Como evento de comparación considere un experimento donde hay canicas en una urna: 5 rojas y 5 blancas. Seleccionamos una canica al azar. Usaremos esta urna como comparación para calibrar nuestra creencia acerca de la tormenta en la autopista. Ahora, considere el siguiente par de apuestas de las cuales puede elegir sólamente una:
A. Obtiene $1000 si hay una tormenta que ocasiona el cierre de la autopista el próximo 20 de noviembre.
B. Obtiene $1000 si seleccionas una canica roja de la urna que contiene 5 canicas rojas y 5 blancas.
Si prefiere la apuesta B quiere decir que considera que la probabilidad de tormenta es menor a 0.5 (50%), por lo que al menos sabe que su creencia subjetiva de una la probabilidad de tormenta es menor a 0.5. Podemos continuar con el proceso para tener una mejor estimación de la creencia subjetiva.
A. Obtiene $1000 si hay una tormenta que ocasiona el cierre de la autopista el próximo 20 de noviembre.
C. Obtiene $1000 si selecciona una canica roja de la urna que contiene 1 canica roja y 9 blancas.
Si ahora selecciona la apuesta A, esto querría decir que considera que la probabilidad de que ocurra una tormenta es mayor a 0.1 (10%). Si consideramos ambas comparaciones tenemos que su probabilidad subjetiva se ubica entre 0.1 y 0.5.
Por tanto, el proceso de calibración nos permite cuantificar nuestras creencias subjetivas. Sin embargo, cuando hay muchos posibles resultados de un evento, es practicamente imposible calibrar las creencias subjetivas para cada resultado; en su lugar, podemos usar una función matemática.
Ejercicio: ¿Cuántos analfabetas dirías que había en México en 2015? Da un intervalo del 90% de confianza para esta cantidad.
Más de calibración:
Prueba de calibración de Messy Matters.
Más pruebas en An Educated Guess.
Ejemplo: podemos pensar que una mujer mexicana promedio mide 156 cm pero estar abierto a la posibilidad de que el promedio sea un poco mayor o menor. En este caso podríamos describir nuestras creencias a través de una curva con forma de campana y centrada en 156.
No olvidemos que estamos describiendo probabilidades, subjetivas o no deben cumplir los axiomas de probabilidad. Es por esto que la función matemática debe conformar una distribuión de probabilidad.
En el caso de parámetros \(\theta\), si \(p(\theta)\) representa el grado de nuestra creencia en los valores de \(\theta\), entonces la media de \(p(\theta)\) se puede pensar como un valor de \(\theta\) que representa nuestra creencia típica o central. Por su parte, la varianza de \(\theta\), que mide que tan dispersa esta la distribución, se puede pensar como la incertidumbre entre los posibles valores.
Thomas Bayes (1702-1761) fue un matemático y ministro de la iglesia presbiteriana, en 1764 se publicó su famoso teorema.
Una aplicación crucial de la regla de Bayes es determinar la probabilidad de un modelo dado un conjunto de datos. Lo que el modelo determina es la probabilidad de los datos condicional a valores particulares de los parámetros y a la estructura del modelo.
Por su parte, podemos utilizar la regla de Bayes para ir de la probabilidad de los datos, dado el modelo, a la probabilidad del modelo, dados los datos.
Ejemplo: Lanzamientos de monedas. Comencemos recordando la regla de Bayes usando dos variables aleatorias discretas. Lanzamos una moneda 3 veces, sea \(X\) la variable aleatoria correspondiente al número de águilas observadas y \(Y\) registra el número de cambios entre águilas y soles.
Escribe la distribución conjunta de las variables, y las distribuciones marginales.
Considera la probabilidad de observar un cambio condicional a que observamos un águila y compara con la probabilidad de observar un águila condicional a que observamos un cambio.
Para entender probabilidad condicional podemos pensar en restringir nuestra atención a una única fila o columna de la tabla. Supongamos que alguien lanza una moneda 3 veces y nos informa que la secuencia contiene exactamente un cambio. Dada esta información podemos restringir nuestra atención a la fila correspondiente a un solo cambio. Sabemos que ocurrió uno de los eventos de esa fila. Las probabilidades relativas de los eventos de esa fila no han cambiado pero sabemos que la probabilidad total debe sumar uno, por lo que simplemente normalizamos dividiendo entre \(p(C=1)\).
Por otro lado, cuando no sabemos nada acerca del número de cambios, todo lo que sabemos del número de águilas está contenido en la distribución marginal de \(X\).
Una de las aplicaciones más importantes de la regla de Bayes es cuando las variables fila \(X\) y columna \(Y\) son datos y parámetros del modelo, respectivamente.
Un modelo especifica la probabilidad de valores particulares de los datos dado la estructura del modelo y valores de los parámetros.
Ejemplo: En un modelo de lanzamientos de monedas tenemos \(p(x = A|\theta)=\theta\) y \(p(x = S|\theta)= 1- \theta\).
De manera general, el modelo especifica: \[p(datos|\text{valores de parámetros y estructura del modelo})\]
y usamos la regla de Bayes para convertir la expresión anterior a lo que nos interesa de verdad, que es, que tanta certidumbre tenemos del modelo condicional a los datos:
\[p(\text{valores de parámetros y estructura del modelo} | datos)\] Una vez que observamos más datos, usamos la regla de Bayes para determinar o actualizar nuestras creencias de posibles parámetros y modelos.
En escencia, el proceso se reduce a calcular:
\[p(\theta|x) = \frac{p(x|\theta)p(\theta)}{p(x)}\]
donde:
Antes de observar datos \(x\), cuantificamos la información (o incertidumbre) acerca del parámetro desconocido \(\theta\) mediante la distribución de probabilidad inicial o a priori \(p(\theta)\). Esto es, la distribución a priori resume nuestras creencias acerca del parámetro ajenas a los datos.
Por otra parte, cuantificamos la información de \(\theta\) asociada a \(x\) mediante la función de verosimilitud \(p(x|\theta) = \mathcal{L}(\theta)\).
Combinamos la información a priori y la información que provee \(x\) mediante el teorema de Bayes obteniendo así la distribución posterior \(p(\theta|x) \propto p(x|\theta)p(\theta)\).
La evidencia \(p(x)\) es la probabilidad de los datos dado el modelo; se determina sumando (caso discreto) o integrando (caso contínuo) \(p(x|\theta)p(\theta)\) a través de todos los valores de los parámetros.
Ejemplo: Ingesta calórica de estudiantes. Supongamos que nos interesa aprender los hábitos alimenticos de los estudiantes universitarios en México, y escuchamos que de acuerdo a investigaciones se recomienda que un adulto promedio ingiera 2500 kcal.
Buscamos conocer qué proporción de los estudiantes siguen esta recomendación. Para ello tomaremos una muestra aleatoria de estudiantes del ITAM.
Denotemos por \(\theta\) la proporción de estudiantes que ingieren en un día 2500 kcal o más. El valor de \(\theta\) es desconocido, y desde el punto de vista bayesiano cuando tenemos incertidumbre de algo (puede ser un parámetro o una predicción) lo vemos como una variable aleatoria y por tanto tiene asociada una distribución de probabilidad que actualizaremos conforme obtenemos información (observamos datos).
theta <- seq(0.05, 0.95, 0.1)
pesos.prior <- c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)
prior <- pesos.prior/sum(pesos.prior)
prior_df <- data_frame(theta, prior = round(prior, 3))
prior_df
#> # A tibble: 10 x 2
#> theta prior
#> <dbl> <dbl>
#> 1 0.0500 0.0350
#> 2 0.150 0.180
#> 3 0.250 0.277
#> 4 0.350 0.249
#> 5 0.450 0.159
#> 6 0.550 0.0730
#> 7 0.650 0.0240
#> 8 0.750 0.00300
#> 9 0.850 0.
#> 10 0.950 0.
\[p(x_1,...,x_n|\theta) = \prod_{i=1}^n p(x_i|\theta)= \theta^z(1 - \theta)^{n-z}\]
donde \(z\) denota el número de estudiantes que ingirió al menos 2500 kcal y \(n-z\) el número de estudiantes que ingirió menos de 2500 kcal.
\[p(\theta|x) = \frac{p(x_1,...,x_N|\theta)p(\theta)}{p(x)} \propto p(\theta)\mathcal{L}(\theta)\]
Continuando con el ejemplo de ingesta calórica de estudiantes universitarios, usamos la inicial discreta que discutimos (tabla de pesos normalizados) y supongamos que tomamos una muestra de 30 alumnos de los cuales \(z=11\) ingieren al menos 2500 kcal. Calculemos la distribución posterior de \(\theta\), con \(0<\theta<1\):
library(LearnBayes)
N <- 30 # estudiantes
z <- 11 # éxitos
# Verosimilitud
Like <- theta ^ z * (1 - theta) ^ (N - z)
product <- Like * prior
# Distribución posterior (normalizamos)
post <- product / sum(product)
dists <- bind_cols(prior_df, post = post)
round(dists, 3)
#> # A tibble: 10 x 3
#> theta prior post
#> <dbl> <dbl> <dbl>
#> 1 0.0500 0.0350 0.
#> 2 0.150 0.180 0.00600
#> 3 0.250 0.277 0.220
#> 4 0.350 0.249 0.529
#> 5 0.450 0.159 0.224
#> 6 0.550 0.0730 0.0210
#> 7 0.650 0.0240 0.
#> 8 0.750 0.00300 0.
#> 9 0.850 0. 0.
#> 10 0.950 0. 0.
# También podemos usar la función pdisc
pdisc(p = theta, prior = prior, data = c(z, N - z)) %>% round(3)
#> [1] 0.000 0.006 0.220 0.529 0.224 0.021 0.000 0.000 0.000 0.000
# Alargamos los datos para graficar
dists_l <- dists %>%
gather(dist, val, prior:post)
ggplot(dists_l, aes(x = theta, y = val)) +
geom_bar(stat = "identity", fill = "darkgray") +
facet_wrap(~ dist) +
labs(x = expression(theta), y = expression(p(theta)))
Ejercicio: ¿Cómo se ve la distribución posterior si tomamos una muestra de tamaño 90 donde observamos la misma proporción de éxitos. Realiza los cálculos y graficala como un panel adicional de la gráfica anterior.
Pregunta: ¿Cómo definirías la distribución inicial si no tuvieras conocimiento de los artículos y expertos?
Ahora, en el teorema de Bayes también encontramos el término \(p(x)\) que denominamos la evidencia, esta última también se conoce como verosimilitud marginal o inicial predictiva.
La evidencia es la probabilidad de los datos de acuerdo al modelo y se calcula sumando o integrando \(p(x, \theta)\) a lo largo de todos los posibles valores de los parámetros \(\theta\).
Es importante notar que hablamos de valores de los parámetros \(\theta\) únicamente en el contexto de un modelo particular \(M\) pues éste es el que da sentido a los parámetros. Podemos hacer evidente el modelo en la notación,
\[p(\theta|x,M)=\frac{p(x|\theta,M)p(\theta|M)}{p(x|M)}\]
En este contexto la evidencia se define como: \[p(x|M)=\sum_{\theta} p(x|\theta,M)p(\theta|M)\] en el caso discreto o, \[p(x|M)=\int p(x|\theta,M)p(\theta|M)d\theta\] en el caso contínuo.
La notación anterior es conveniente sobre todo cuando estamos considerando más de un modelo y queremos usar los datos para determinar la certeza que tenemos en cada modelo.
Supongamos que tenemos dos modelos \(M_1\) y \(M_2\), entonces podemos calcular el cociente de \(p(M_1|x)\) y \(p(M_2|x)\) obteniendo:
\[\frac{p(M_1|x)}{p(M_2|x)} = \frac{p(x|M_1) \cdot p(M_1)}{p(x|M_2)\cdot p(M_2)}\]
El cociente de evidencias \(\frac{p(x|M_1)}{p(x|M_2)}\) se conoce como factor de Bayes.
Vimos que la regla de Bayes nos permite pasar del conocimiento inicial \(p(\theta)\) al posterior \(p(\theta|x)\) conforme recopilamos datos. Supongamos ahora que observamos más datos, los denotamos \(x'\), podemos volver a actualizar nuestras creencias pasando de \(p(\theta|x)\) a \(p(\theta|x,x')\). Entonces podemos preguntarnos si nuestro conocimiento posterior cambia si actualizamos de acuerdo a \(x\) primero y \(x'\) después o vice-versa. La respuesta es que si \(p(x|\theta)\) y \(p(x'|\theta)\) son iid entonces el orden en que actualizamos nuestro conocimiento no afecta la distribución posterior. Esta propiedad se denomina invarianza al orden.
La invarianza al orden tiene sentido intuitivamente: Si la función de verosimilitud no depende del tiempo o del ordenamineto de los datos, entonces la posterior tampoco tiene porque depender del ordenamiento de los datos.
Los tres objetivos de la inferencia son: estimación de parámetros, predicción de valores y comparación de modelos.
Ejemplo: Considere que proponemos los siguientes dos modelos (distribución inicial y función de verosimilitud) para una muestra que consta de una susesión de \(N\) lanzamientos de moneda en la que salieron \(z\) águilas:
# Modelo 1, creamos una inicial que puede tomar 3 valores
p_M1 <- data_frame(theta = c(0.25, 0.5, 0.75), prior = c(0.25, 0.5, 0.25),
modelo ="M1")
p_M1
#> # A tibble: 3 x 3
#> theta prior modelo
#> <dbl> <dbl> <chr>
#> 1 0.250 0.250 M1
#> 2 0.500 0.500 M1
#> 3 0.750 0.250 M1
# Modelo 2, Creamos una inicial que puede tomar 51 valores
p <- seq(0, 24, 1)
p
#> [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
#> [24] 23 24
p2 <- c(p, 25, sort(p, decreasing = TRUE))
p2
#> [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
#> [24] 23 24 25 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5
#> [47] 4 3 2 1 0
p_M2 <- data_frame(theta = seq(0, 1, 0.02), prior = p2 / sum(p2),
modelo = "M2")
p_M2
#> # A tibble: 51 x 3
#> theta prior modelo
#> <dbl> <dbl> <chr>
#> 1 0. 0. M2
#> 2 0.0200 0.00160 M2
#> 3 0.0400 0.00320 M2
#> 4 0.0600 0.00480 M2
#> 5 0.0800 0.00640 M2
#> 6 0.100 0.00800 M2
#> 7 0.120 0.00960 M2
#> 8 0.140 0.0112 M2
#> 9 0.160 0.0128 M2
#> 10 0.180 0.0144 M2
#> # ... with 41 more rows
#Parámetros
N <- 20 # estudiantes
z <- 12 # éxitos
dists_h <- bind_rows(p_M1, p_M2) %>% # base de datos horizontal
group_by(modelo) %>%
mutate(
Like = theta ^ z * (1 - theta) ^ (N - z), # verosimilitud
posterior = (Like * prior) / sum(Like * prior)
)
dists_h
#> # A tibble: 54 x 5
#> # Groups: modelo [2]
#> theta prior modelo Like posterior
#> <dbl> <dbl> <chr> <dbl> <dbl>
#> 1 0.250 0.250 M1 5.97e- 9 2.49e- 3
#> 2 0.500 0.500 M1 9.54e- 7 7.96e- 1
#> 3 0.750 0.250 M1 4.83e- 7 2.02e- 1
#> 4 0. 0. M2 0. 0.
#> 5 0.0200 0.00160 M2 3.48e-21 9.53e-18
#> 6 0.0400 0.00320 M2 1.21e-17 6.62e-14
#> 7 0.0600 0.00480 M2 1.33e-15 1.09e-11
#> 8 0.0800 0.00640 M2 3.53e-14 3.86e-10
#> 9 0.100 0.00800 M2 4.30e-13 5.89e- 9
#> 10 0.120 0.00960 M2 3.21e-12 5.26e- 8
#> # ... with 44 more rows
dists <- dists_h %>% # base de datos larga
gather(dist, valor, prior, Like, posterior) %>%
mutate(dist = factor(dist, levels = c("prior", "Like", "posterior")))
dists
#> # A tibble: 162 x 4
#> # Groups: modelo [2]
#> theta modelo dist valor
#> <dbl> <chr> <fct> <dbl>
#> 1 0.250 M1 prior 0.250
#> 2 0.500 M1 prior 0.500
#> 3 0.750 M1 prior 0.250
#> 4 0. M2 prior 0.
#> 5 0.0200 M2 prior 0.00160
#> 6 0.0400 M2 prior 0.00320
#> 7 0.0600 M2 prior 0.00480
#> 8 0.0800 M2 prior 0.00640
#> 9 0.100 M2 prior 0.00800
#> 10 0.120 M2 prior 0.00960
#> # ... with 152 more rows
ggplot(filter(dists, modelo == "M1"), aes(x = theta, y = valor)) +
geom_bar(stat = "identity", fill = "darkgray") +
facet_wrap(~ dist, scales = "free") +
scale_x_continuous(expression(theta), breaks = seq(0, 1, 0.2)) +
labs(y = "")
ggplot(filter(dists, modelo == "M2"), aes(x = theta, y = valor)) +
geom_bar(stat = "identity", fill = "darkgray") +
facet_wrap(~ dist, scales = "free") +
scale_x_continuous(expression(theta), breaks = seq(0, 1, 0.2)) +
labs(y = "")
Si nuestro conocimiento actual consta únicamente de las creencias iniciales de la distribución inical \(p(\theta)\), utilizamos la distribución predictiva inicial: \[p(y) =\sum_{\theta} p(y|\theta)p(\theta)\] en el caso discreto, o \[p(y) =\int p(y|\theta)p(\theta)d\theta\] en el caso contínuo.
Notese que la ecuación anterior coincide con la correspondiente a la evidencia, con la diferencia de que la evidencia se refiere a un valor observado y en esta ecuación estamos calculando la probabilidad de cualquier valor \(y\).
Si nuestro conocimiento actual involucra datos observados \(x\), utilizamos la distribución predictiva posterior: \[p(y|x) =\sum_{\theta} p(y|\theta)p(\theta|x)\] en el caso discreto, o \[p(y|x) =\int p(y|\theta)p(\theta|x)d\theta\] en el caso contínuo.
Por ejemplo podemos usar las creencias iniciales del modelo 1, que propusimos arriba para calcular la probabilidad predictiva de observar águila:
\[p(y=A) = \sum_{\theta}p(y=A|\theta)p(\theta) = 0.5\] Vale la pena destacar que las prediciones son probabilidades de cada posible valor condicional a nuestro modelo de creencias actuales. Si nos interesa predecir un valor particular en lugar de de una distribución a lo largo de todos los posibles valores podemos usar la media de la distribución predictiva. Por tanto el valor a predecir sería (en el caso contínuo):
\[E[y|x]=\int y \ p(y|x) dy\]
La integral anterior únicamente tiene sentido si \(y\) es una variable continua. Si \(y\) es nominal, entonces podemos usar el valor más probable.
Ejercicio: Calcula las probabilidades predictivas para \(y=0,1\) de los modelos 1 y 2 del ejemplo anterior usando la distribución posterior de cada modelo.
Ejemplo: El primero modelo del ejemplo anterior supusimos que el parámetro \(\theta\) únicamente puede tomar uno de 3 valores (0.25, 0.5, 0.75); esta restricción dió lugar a un modelo simple. Por su parte, el modelo 2 es más complejo y permite muchos más valores de \(\theta\) (51). La forma de la distribución inicial es triangular en ambos casos, el valor de mayor probabilidad inicial es \(\theta = 0.50\) y reflejamos que creemos que es menos factible que el valor se encuentre en los extremos.
Podemos calcular el factor de Bayes para distintos datos observados:
factorBayes <- function(N, z){
evidencia <- bind_rows(p_M1, p_M2) %>% # base de datos horizontal
group_by(modelo) %>%
mutate(
Like = theta ^ z * (1 - theta) ^ (N - z), # verosimilitud
posterior = (Like * prior) / sum(Like * prior)
) %>%
summarise(evidencia = sum(prior * Like))
print(evidencia)
return(evidencia[1, 2] / evidencia[2, 2])
}
factorBayes(50, 25)
#> # A tibble: 2 x 2
#> modelo evidencia
#> <chr> <dbl>
#> 1 M1 4.44e-16
#> 2 M2 2.76e-16
#> evidencia
#> 1 1.61
factorBayes(100, 75)
#> # A tibble: 2 x 2
#> modelo evidencia
#> <chr> <dbl>
#> 1 M1 9.46e-26
#> 2 M2 4.16e-26
#> evidencia
#> 1 2.27
factorBayes(100, 10)
#> # A tibble: 2 x 2
#> modelo evidencia
#> <chr> <dbl>
#> 1 M1 1.36e-18
#> 2 M2 2.47e-16
#> evidencia
#> 1 0.0055
factorBayes(40, 38)
#> # A tibble: 2 x 2
#> modelo evidencia
#> <chr> <dbl>
#> 1 M1 0.000000279
#> 2 M2 0.00000894
#> evidencia
#> 1 0.0313
¿Cómo explicarías los resultados anteriores? En el modelo complejo \(\theta\) puede tomar más valores, entonces si en una sucesión de volados observamos 10% de águilas, el modelo simple no cuenta con un valor de \(\theta\) cercano al resultado observado, pero el modelo complejo si. Por otra parte, para valores de \(\theta\) que se encuentran en ambos modelos la probabilidad inicial de esos valores es mayor en el caso del modelo simple. Por lo tanto si los datos observados resultan en valores de \(\theta\) congruentes con el modelo simple, creeríamos en el modelo simple más que en el modelo más complicado.
Es importante destacar que la comparación de modelos nos habla únicamente de la evidencia relativa de dos modelos; sin embargo, puede que ninguno de los modelos que estamos considerando sean adecuados para nuestros datos, por lo que más adelante estudiaremos maneras de evaluar un modelo.
En la inferencia Bayesiana se requiere calcular el denominador de la regla de Bayes \(p(x)\); es común que esto requiera que se calcule una integral complicada, que en ocasiones no tenga una forma cerrada. Sin embargo, hay algunas maneras de evitar esto:
El camino tradicional consiste en usar funciones de verosimilitud con distribuciones iniciales conjugadas. Cuando una distribución inicial es conjugada de la función de verosimilitud, la integral resulta en una distribución posterior con la misma forma funcional que la distribución inicial.
Otra alternativa es aproximar la integral numéricamente. Cuando el espacio de parámetros es de dimensión chica, se puede aproximar con una cuadrícula de puntos y la integral se puede calcular sumando a través de dicha cuadrícula. Sin embargo, cuando el espacio de parámetros aumenta de dimensión el número de puntos necesarios para la aproximación crece demasiado y hay que recurrir a otas técnicas.
Se han desarrollado una clase de métodos de simulación para poder calcular la distribución posterior, estos se conocen como cadenas de Markov via Monte Carlo (MCMC por sus siglas en inglés). El desarrollo de los métodos MCMC es lo que ha propiciado el desarrollo de la estadística bayesiana en años recientes.
Para ilustrar estas 3 aproximaciones a la integral que define la distribución posterior, consideraremos el ejemplo de una distribución Bernoulli.
Comenzaremos con el modelo Beta-Binomial. Recordemos que si \(X\) es una variable aleatoria con dos posibles resultados, entonces \(X\) se distribuye Bernoulli y la función de probabilidad esta definida por:
\[p(x|\theta)=\theta^x(1-\theta)^{1-x}\]
Usualmente pensamos en la fórmula anterior como una función de \(x\), donde \(x\) toma uno de dos valores: 0 o 1. También podemos pensar que \(x\) esta fija (observada) y la expresión es una función de \(\theta\), en este caso, la ecuación especifica la probabilidad de un valor fijo \(x\) para un valor de \(\theta\) y la llamamos la función de verosimilitud de \(\theta\). En inferencia Bayesiana, la función \(p(x|\theta)\) usualmente se considera como función de \(\theta\).
Ahora, si lanzamos una moneda \(N\) veces tenemos un conjunto de datos \(\{x_1,...,x_N\}\). Si suponemos que los lanzamientos son i.i.d., la probabilidad de observar el conjunto de \(N\) lanzamientos es el producto de las probabilidades de cada observación:
\[p(x_1,...,x_N|\theta) = \prod_{n=1}^N p(x_n|\theta)\] \[= \theta^z(1 - \theta)^{N-z}\] donde \(z\) denota el número de éxitos (águilas).
Ahora, en principio, para describir nuestras creencias iniciales podríamos usar cualquier función de densidad de probabilidad inicial \(p(\theta)\) con soporte en \([0, 1]\). Sin embargo, sería conveniente que el producto \(p(x|\theta)p(\theta)\) (el numerador de la fórmula de Bayes) resulte en una función con la misma forma funcional que \(p(\theta)\). Cuando éste es el caso, las distribuciones inicial y posterior se describen con la misma fisma familia de distribuciones. Esto es conveninte pues si obtenemos nueva información podemos actualizar nuestro conocimiento de manera inmediata, conservando la forma de las distribuciones.
Cuando las funciones \(p(x|\theta)\) y \(p(\theta)\) se combinan de tal manera que la distribución posterior pertenece a la misma familia (tiene la misma forma funcional) que la distribución inicial, entonces decimos que \(p(\theta)\) es conjugada con \(p(x|\theta)\).
Vale la pena notar que la distribución inicial es conjugada con una función de verosimilitud particular.
Una distribución conjugada con \(p(x|\theta) = \theta^z(1 - \theta)^{N-z}\) es una \(Beta(a, b)\) \[p(\theta) = \frac {\theta^{a-1}(1-\theta)^{b-1}}{B(a,b)}\]
Para describir nuestro conocimiento inicial podemos explorar la media y la varianza de la distribución beta. La media es \[\bar{\theta} = a/(a+b)\] por lo que si \(a=b\) la media es 0.5 y conforme aumenta \(a\) en relación a \(b\) aumenta la media. La varianza es:
\[s^2 = \frac{\bar{\theta}(1-\bar{\theta})}{(a+b+1)}\]
Una manera de seleccionar los parámetros \(a\) y \(b\) es pensar en la proporción media de águilas (m) y el tamaño de la muestra (n). Ahora, \(m=a/(a+b)\) y \(n = a+b\), por lo que obtenemos
\[a=mn, \quad b=(1-m)n\]
Otra manera es comenzar con la media y la varianza. Al usar este enfoque debemos recordar que la desviación estándar \(s\) debe tener sentido en el contexto de la densidad beta. En particular la desviación estándar típicamente es menor a 0.289, la cual corresponde a la desviación estándar de una distribución uniforme contínua. Entonces, para una distribución beta con media \(\bar{\theta}\) y varianza \(s^2\), los parámetros son: \[a=\bar{\theta}\bigg(\frac{\bar{\theta}(1-\bar{\theta})}{s^2}- 1\bigg), \quad b=(1-\bar{\theta})\bigg(\frac{\bar{\theta}(1-\bar{\theta})}{s^2}- 1\bigg)\]
Una vez que hemos determinado una distribución inicial conveniente para la función de verosimilitud Bernoulli, calculamos la distribución posterior.
Supongamos que observamos \(N\) lanzamientos de los cuales \(z\) son águilas, entonces podemos ver que la distribución posterior es nuevamente una densidad Beta.
\[p(\theta|z)\propto \frac{\theta^{a+z-1}(1 -\theta)^{(N-z+b)-1}}{B(a,b)} = p(x|\theta) p(\theta)\]
Concluímos entonces que si la distribución inicial es \(Beta(a,b)\), la posterior es \(Beta(z+a,N-z+b)\).
Vale la pena explorar la relación entre las medias de la distribución inicial y la distribución posterior. La media incial es \[\frac{a}{(a+b)}\] y la media posterior es \[\frac{(z+a)}{(z+a) + (N-z+b)}=\frac{z+a}{N+a+b}\] Podemos hacer algunas manipulaciones algebráicas para escribirla como:
\[\frac{z+a}{N+a+b}=\frac{z}{N}\frac{N}{N+a+b} + \frac{a}{a+b}\frac{a+b}{N+a+b}\]
es decir, podemos escribir la media posterior como un promedio ponderado entre la media inicial \(a/(a+b)\) y la proporción observada \(z/N\).
Ejemplo: Consideremos dos distribuciones iniciales beta diferentes:
N = 14; z = 11; a = 1; b = 1
base <- ggplot(data_frame(x = c(0, 1)), aes(x))
p1 <- base +
stat_function(fun = dbeta, args = list(shape1 = a, shape2 = b),
aes(colour = "inicial"), show.legend = FALSE) +
stat_function(fun = dbeta, args = list(shape1 = z + 1, shape2 = N - z + 1),
aes(colour = "verosimilitud"), show.legend = FALSE) +
stat_function(fun = dbeta, args = list(shape1 = a + z, shape2 = N - z + b),
aes(colour = "posterior"), show.legend = FALSE) +
labs(y = "", colour = "", x = expression(theta))
N = 14; z = 11; a = 100; b = 100
p2 <- base +
stat_function(fun = dbeta, args = list(shape1 = a, shape2 = b),
aes(colour = "inicial")) +
stat_function(fun = dbeta, args = list(shape1 = z + 1, shape2 = N - z + 1),
aes(colour = "verosimilitud")) +
stat_function(fun = dbeta, args = list(shape1 = a + z, shape2 = N - z + b),
aes(colour = "posterior")) +
labs(y = "", colour = "", x = expression(theta))
grid.arrange(p1, p2, nrow = 1, widths = c(0.38, 0.62))
Una manera de resumir la distribución posterior es a través de intervalos de probabilidad. En inferencia Bayesiana, los intervalos de probabilidad se utilizan para establecer qué valores del parámetro son creíbles, por lo que se denominan intervalos de credibilidad.
Ejercicio: Calcula un intervalo de credibilidad del 95% para cada una de las distribuciones posteriores del ejemplo anterior.
\[p(y = 1|x) = \int p(y=1|\theta)p(\theta|x)d\theta\] \[=\int \theta \ p(\theta|z,N) d\theta\] \[=\frac{z+a}{N+a+b}\] Esto es, la probabilidad predictiva de águila es la media de la distribución posterior sobre \(\theta\).
\[p(x|M)=\int p(x|\theta,M)p(\theta|M)d\theta\]
En este caso, los datos se resumen mediante las cantidades \(z\) y \(N\). Para una distribución incial beta es fácil calcular la evidencia:
\[p(x|M) \equiv p(z,N)=\frac{B(z+a,N-z+b)}{B(a,b)}\]
En nuestro ejemplo:
# N = 14, z = 11, a = 1, b = 1
P_M1 <- beta(12, 4) / beta(1, 1)
P_M1
#> [1] 0.000183
# N = 14, z = 12, a = 100, b = 100
P_M2 <- beta(126, 126) / beta(100, 100)
P_M2
#> [1] 1.98e-16
factor_Bayes <- P_M1 / P_M2
factor_Bayes
#> [1] 9.26e+11
Vemos que el modelo 1 se ajusta mucho mejor a los datos: la distribución inicial del modelo 2 tiene un pico en 0.5 mientras que la otra es uniforme. Debido a que la proporción de águilas observadas en la muestra no es cercana a 0.5, la distribución inicial picuda (modelo 2) no captura bien los datos.
En general, preferimos un modelo con un valor mayor de \(p(x|\theta)\), pero la preferencia no es absoluta ya que una diferencia chica no nos dice mucho; debemos considerar que los datos no son mas que una muestra aleatoria.
Ejercicio: Supongamos que nos interesa analizar el IQ de una muestra de estudiantes del ITAM y suponemos que el IQ de un estudiante tiene una distribución normal \(x \sim N(\theta, \sigma^2)\) con \(\sigma ^ 2\) conocida.
Considera que observamos el IQ de un estudiante \(x\). La función de verosimilitud del modelo es: \[p(x|\theta)=\frac{1}{\sqrt{2\pi\sigma^2}}exp\left(-\frac{1}{2\sigma^2}(x-\theta)^2\right)\]
Realizaremos un análisis bayesiano, por lo que hace falta establecer una distribución inicial. Elegiremos que \(p(\theta)\) se distribuya \(N(\mu, \tau^2)\) donde elegimos los parámetros \(\mu, \tau\) que mejor describan nuestras creencias iniciales.
Calcula la distribución posterior \(p(\theta|x) \propto p(x|\theta)p(\theta)\), utilizando la distribución inicial y la función de verosimilitud que definimos arriba. Una vez que realices la multiplicación debes identificar el núcleo de una distribución Normal, ¿cuáles son los parámetros (media y varianza) de la distribución posterior?
Supongamos que la distribución beta no describe nuestras creencias de manera adecuada. Por ejemplo, mis creencias podrían estar mejor representadas por una distribución trimodal: con tres máximos (ej. la moneda está fuertemente sesgada hacia sol, fuertemente sesgada hacia águila o es justa). No hay parámetros en una beta que puedan describir este patrón, y probablemente no esposible realizar la integración analítica para encontrar la distribución posterior con una función que describe una distribución inicial de este tipo.
Exploraremos entonces una técnica de aproximación numérica de la distribución posterior que consiste en definir la distribución inicial en una cuadrícula de valores de \(\theta\). En este método no necesitamos describir nuestras creencias mediante una función matemática ni realizar integración analítica.
Suponemos que existe únicamente un número finito de valores de \(\theta\) que creemos que pueden ocurrir (en el primer ejemplo que estudiamos usamos esta técnica). Es así que la regla de Bayes se escribe como:
\[p(x|\theta)=\frac{p(x|\theta)p(\theta)}{\sum_{\theta}p(x|\theta)p(\theta)}\]
Esta técnica nos permite aproximar la integral que aparecería en el denominador anterior discretizando una distribución inicial continua mediante una cuadrícula de masas de probabilidad discreta y por tanto podemos usar la versión discreta de la regla de Bayes. El proceso consiste en dividir el dominio en regiones, crear un rectángulo con la altura correspondiente al valor de la densidad en el punto medio. Aproximamos el área de cada región mediante la altura del rectángulo.
Ejemplo: Distribución inicial contínua uniforme. Si discretizamos la distribución considerando 20 valores del parámetro \(\theta\) obtenemos:
# N = 14, z = 11, a = 1, b = 1
N = 14; z = 11
inicial <- data.frame(theta = seq(0.05, 1, 0.05), inicial = rep(1/20, 20))
dists_h <- inicial %>%
mutate(
verosimilitud = theta ^ z * (1 - theta) ^ (N - z), # verosimilitud
posterior = (verosimilitud * inicial) / sum(verosimilitud * inicial)
)
dists <- dists_h %>% # base de datos larga
gather(dist, valor, inicial, verosimilitud, posterior) %>%
mutate(dist = factor(dist, levels = c("inicial", "verosimilitud", "posterior")))
ggplot(dists, aes(x = theta, y = valor)) +
geom_point() +
facet_wrap(~ dist, scales = "free") +
scale_x_continuous(expression(theta), breaks = seq(0, 1, 0.2)) +
labs(y = "")
y lo podemos comparar con la versión continua (distribución beta con \(a=b=1\)):
Podemos calcular resúmenes estadísticos de la distribución de \(\theta\) de manera aproximada. Por ejemplo, la media se aproxima mediante promedio ponderado por las probabilidades:
\(\bar{\theta}=\sum_{\theta} \theta p(\theta|x)\)
head(dists_h)
#> theta inicial verosimilitud posterior
#> 1 0.05 0.05 4.19e-15 1.14e-12
#> 2 0.10 0.05 7.29e-12 1.99e-09
#> 3 0.15 0.05 5.31e-10 1.45e-07
#> 4 0.20 0.05 1.05e-08 2.86e-06
#> 5 0.25 0.05 1.01e-07 2.75e-05
#> 6 0.30 0.05 6.08e-07 1.66e-04
sum(dists_h$posterior * dists_h$theta)
#> [1] 0.75
En cuanto a los intervalos de probabilidad, debido a que estamos usando masas discretas, la suma de las masas en un intervalo usualmente no será exactamente igual al coeficiente de confianza. Para un coeficiente de confianza de 95%, tenemos que elegir los puntos tales que la masa sea mayor a igual a 95% y la masa total sea lo menor posible. En nuestro ejemplo podemos usar cuantiles:
dist_cum <- cumsum(dists_h$posterior) #vector de distribución acumulada
lb <- which.min(dist_cum < 0.05) - 1
ub <- which.min(dist_cum < 0.975)
dists_h$theta[lb]
#> [1] 0.5
dists_h$theta[ub]
#> [1] 0.9
\[p(y|x)=\int p(y|\theta)p(\theta|x)d\theta\] \[\approx \sum_{\theta} p(y|\theta)p(\theta|x)\]
Ejercicio: Calcula la probabilidad predictiva para \(y=1\) usando los datos del ejemplo anterior.
\[p(x|M)=\int p(x|\theta,M)p(\theta|M)d\theta\]
se convierte en una suma
\[p(x|M)\approx \sum_{\theta} p(x|\theta,M)p(\theta|M)d\theta\]
Para el ejemplo del experimento Bernoulli donde definimos los modelos \(M_1\) y \(M_2\) tendríamos:
# calcula el factor de Bayes para el experimento Bernoulli Modelos M1 y M2
factorBayes <- function(M, s){
evidencia <- rbind(p_M1, p_M2) %>% # base de datos horizontal
group_by(modelo) %>%
mutate(
Like = theta ^ s * (1 - theta) ^ (M - s), # verosimilitud
posterior = (Like * prior) / sum(Like * prior)
) %>%
summarise(evidencia = sum(prior * Like))
print(evidencia)
return(evidencia[1, 2] / evidencia[2, 2])
}
factorBayes(50, 25)
#> # A tibble: 2 x 2
#> modelo evidencia
#> <chr> <dbl>
#> 1 M1 4.44e-16
#> 2 M2 2.76e-16
#> evidencia
#> 1 1.61
Hay ocasiones en las que los métodos de inicial conjugada y aproximación por cuadrícula no funcionan; hay casos en los que la distribución conjugada no describe nuestras creencias iniciales o la aproximación numérica no es factible debido al gran número de parámetros involucrados. Es por ello que surge la necesidad de utilizar métodos de Monte Carlo vía Cadenas de Markov (MCMC).
El algoritmo de Metropolis asume que podemos calcular \(p(\theta)\) para un valor particular de \(\theta\) y el valor de la verosimilitud \(p(x|\theta)\) para cualquier \(x\), \(\theta\) dados. En realidad, el método únicamente requiere que se pueda calcular el producto de la distribución inicial y la función de verosimilitud. Lo que el método produce es una aproximación de la distribución posterior \(p(\theta|x)\) mediante una muestra de valores de \(\theta\) obtenido de dicha distribución.
Ejemplo: Caminata aleatoria. Con el fin de entender el algoritmo comenzaremos estudiando el concepto de caminata aleatoria a través de un ejemplo.
Supongamos que un vendedor de enciclopedias trabaja a lo largo de una cadena finita de islas. Constantemente viaja entre las islas ofreciendo sus productos. Al final de un día de trabajo decide si permanece en la misma isla o se transporta a una de las 2 islas vecinas. El vendedor desconoce a priori la distribución de la población en las islas; sin embargo, una vez que se encuentra en una isla puede investigar la población de la misma y también de la isla a la que se propone viajar después.
El objetivo del vendedor es visitar las islas de manera proporcional a la población de cada una. Con esto en mente el vendedor utiliza el siguiente proceso:
Lanza un volado, si el resultado es águila se propone ir a la isla del lado izquierdo de su ubicación actual y si es sol a la del lado derecho.
Si la isla propuesta en el paso anterior tiene población mayor a la población de la isla actual, el vendedor decide viajar a ella. Si la isla tiene población menor, entonces visita la isla propuesta con una probabilidad que depende de la población relativa de la isla propuesta y la isla actual.
Sea \(\theta\) el número de isla. Sea \(P(\theta_{prop})\) la población de la isla propuesta y \(P(\theta_{actual})\) la población de la isla actual. Entonces el vendedor cambia de isla con probabilidad \(p_{mover}=P(\theta_{prop})/P(\theta_{actual})\).
A la larga, si el vendedor sigue la heurística anterior la probabilidad de que el vendedor se encuentre en alguna de las islas \(P(\theta)\) coincide con la población relativa (normalizada) de la isla.
Consideremos el caso en el que se tienen 10 islas y la diferencia de población entre islas consecutivas es constante:
islas <- data.frame(islas = 1:10, pob = 1:10)
caminaIsla <- function(i){ # i: isla actual
u <- runif(1) # volado
v <- ifelse(u < 0.5, i - 1, i + 1) # isla vecina (índice)
if(v < 1 | v > 10){ # si estas en los extremos y el volado indica salir
return(i)
}
u2 <- runif(1)
p_move = min(islas$pob[v] / islas$pob[i], 1)
if(p_move > u2){
return(v) # isla destino
}
else{
return(i) # me quedo en la misma isla
}
}
pasos <- 100000
camino <- numeric(pasos)
camino[1] <- sample(1:10, 1) # isla inicial
for(j in 2:pasos){
camino[j] <- caminaIsla(camino[j - 1])
}
caminata <- data.frame(pasos = 1:pasos, isla = camino)
plot_caminata <- ggplot(caminata[1: 1000, ], aes(x = pasos, y = isla)) +
geom_point(size = 0.8) +
geom_path(alpha = 0.5) +
coord_flip() +
labs(title = "Caminata aleatoria") +
scale_y_continuous(expression(theta), breaks = 1:10) +
scale_x_continuous("Tiempo")
plot_dist <- ggplot(caminata, aes(x = isla)) +
geom_histogram() +
scale_x_continuous(expression(theta), breaks = 1:10) +
labs(title = "Distribución objetivo (no normalizada)",
y = expression(P(theta)))
grid.arrange(plot_caminata, plot_dist, ncol = 1, heights = c(4, 2))
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
La distribución que deseamos conocer se conoce como distribución objetivo \(P(\theta)\). En este ejemplo sencillo la distribución objetivo es la propabilidad de que el vendedor se encuentre en una isla, la cual por construcción coincide con la población relativa de cada isla (normalizada). Para obtener una buena aproximación de la distribución objetivo debemos permitir que el vendedor recorra las islas durante una sucesión larga de pasos y registramos sus visitas (frecuencias relativas).
Nuestra aproximación de la distribución objetivo es justamente el registro de sus visitas (normalizado). Más aún, debemos tener cuidado y excluir la porción de las visitas que se encuentran bajo la influencia de la posición inicial. Esto es, debemos excluir el periodo de calentamiento.
Una vez que tenemos un registro largo de los viajes del vendedor (excluyendo el calentamiento) podemos aproximar la distribución objetivo \(P(\theta)\) de cada valor de \(\theta\) (el número de isla) simplemente contando el número relativo de veces que el vendedor visitó dicha isla y normalizando.
t <- c(1:10, 20, 50, 100, 200, 1000, 5000)
plots_list <- lapply(t, function(i){
ggplot(caminata[1:i, ], aes(x = isla)) +
geom_histogram() +
labs(y = "", x = "", title = paste("t = ", i, sep = "")) +
scale_x_continuous(expression(theta), breaks = 1:10, limits = c(0, 11))
})
args.list <- c(plots_list,list(nrow=4,ncol=4))
do.call(grid.arrange, args.list)
Escribamos el algoritmo, para esto indexamos las islas por el valor \(\theta\), es así que la isla del extremo izquierdo corresponde a \(\theta=1\) y la población relativa normalizada de cada isla es \(P(\theta)\):
El rango de los posibles valores para moverse, y la probabilidad de proponer cada uno se conoce como distribución propuesta; en nuestro ejemplo sólo toma dos valores cada uno con probabilidad 0.5.
\[p_{mover}=min\bigg( \frac{P(\theta_{propuesta})}{P(\theta_{actual})},1\bigg)\]
Notemos que la distribución objetivo \(P(\theta)\) no necesita estar normalizada, esto es porque lo que nos interesa es el cociente \(P(\theta_{propuesta})/P(\theta_{actual})\).
Entonces, para utilizar el algoritmo necesitamos ser capaces de:
Generar un valor de la distribución propuesta (para crear \(\theta_{propuesta}\)).
Evaluar la distribución objetivo en cualquier valor propuesto (para calcular \(P(\theta_{propuesta})/P(\theta_{actual})\)).
Generar un valor uniforme (para movernos con probabilidad \(p_{mover}\))
Las 3 puntos anteriores nos permiten generar muestras aleatorias de la distribución objetivo, sin importar si esta está normalizada.
Esta técnica es particularmente útil cuando la distribución objetivo es la distribución posterior \(p(\theta|x)\), proporcional a \(p(x|\theta)p(\theta)\).
Para entender porque funciona el algoritmo de Metrópolis hace falta entender 2 puntos:
library(expm)
transMat <- function(P){ # recibe vector de probabilidades (o población)
T <- matrix(0, 10, 10)
n <- length(P - 1) # número de estados
for(j in 2:n - 1){ # llenamos por fila
T[j, j - 1] <- 0.5 * min(P[j - 1] / P[j], 1)
T[j, j] <- 0.5 * (1 - min(P[j - 1] / P[j], 1)) +
0.5 * (1 - min(P[j + 1] / P[j], 1))
T[j, j + 1] <- 0.5 * min(P[j + 1] / P[j], 1)
}
# faltan los casos j = 1 y j = n
T[1, 1] <- 0.5 + 0.5 * (1 - min(P[2] / P[1], 1))
T[1, 2] <- 0.5 * min(P[2] / P[1], 1)
T[n, n] <- 0.5 + 0.5 * (1 - min(P[n - 1] / P[n], 1))
T[n, n - 1] <- 0.5 * min(P[n - 1] / P[n], 1)
T
}
T <- transMat(islas$pob)
w <- c(0, 1, rep(0, 8))
t <- c(1:10, 20, 50, 100, 200, 1000, 5000)
expT <- map_df(t, ~data.frame(t = ., w %*% (T %^% .)))
expT_long <- expT %>%
gather(theta, P, -t) %>%
mutate(theta = parse_number(theta))
ggplot(expT_long, aes(x = theta, y = P)) +
geom_bar(stat = "identity", fill = "darkgray") +
facet_wrap(~ t) +
scale_x_continuous(expression(theta), breaks = 1:10, limits = c(0, 11))
inicioP <- function(i){
w <- rep(0, 10)
w[i] <- 1
t <- c(1, 10, 50, 100)
expT <- map_df(t, ~data.frame(t = ., inicio = i, w %*% (T %^% .))) %>%
gather(theta, P, -t, -inicio) %>%
mutate(theta = parse_number(theta))
expT
}
expT <- map_df(c(1, 3, 5, 9), inicioP)
ggplot(expT, aes(x = as.numeric(theta), y = P)) +
geom_bar(stat = "identity", fill = "darkgray") +
facet_grid(inicio ~ t) +
scale_x_continuous(expression(theta), breaks = 1:10, limits = c(0, 11))
Es importante recalcar que el algoritmo converge a la distribución objetivo sin importar si el punto de partida es una sola isla o cualquier distribución de probabilidad de iniciar en diferentes islas.
En la sección anterior implementamos el algoritmo de Metrópolis en un caso sencillo: las posiciones eran discretas, en una dimensión y la distribución propuesta involucraba únicamente mover a la izquierda o a la derecha. El algoritmo general aplica para valores continuos, en cualquier número de dimensiones y con distribuciones propuesta más generales. Lo esencial del método no cambia para el caso general, esto es:
Tenemos una distribución objetivo \(P(\theta)\) de la cual buscamos generar muestras. Debemos ser capaces de calcular el valor de \(P(\theta)\) para cualquier valor candidato \(\theta\). La distribución objetivo \(P(\theta)\) no tiene que estar normalizada; típicamente \(P(\theta)\) es la distribución posterior de \(\theta\) no normalizada, es decir, es el producto de la verosimilitud y la inicial.
La muestra de la distribución objetivo se genera mediante una caminata aleatoria a través del espacio de parámetros. La caminata inicia en un lugar arbitrario (definido por el usuario). El punto inicial debe ser tal que \(P(\theta)>0\). La caminata avanza en cada tiempo proponiendo un movimiento a una nueva posición mediante una distribución propuesta (1 elemento \(\theta\) muestreado) y después decidiendo si se acepta o no el valor propuesto. Las distribución propuesta puede tener muchas formas, el objetivo es que la distribución propuesta explore el espacio de parámetros de manera eficiente.
Una vez que tenemos un valor propuesto decidimos si aceptar el movimiento calculando:
\[p_{mover}=min\bigg( \frac{P(\theta_{propuesta})}{P(\theta_{actual})},1\bigg)\]
Si el movimiento no es aceptado, el valor \(\theta_{actual}\) se repite en la muestra mientras que el valor \(\theta_{propuesta}\) es descartado de la muestra. Al final obtenemos valores representativos de la distribución objetivo \(\{\theta_1,...,\theta_n\}\).
Es importante recordar que debemos excluir las primeras observaciones pues estas siguen bajo la influencia del valor inicial.
Retomemos el problema de inferencia Bayesiana y veamos como usar el algoritmo de Metrópolis cuando la distribución objetivo es la distribución posterior.
Ejemplo: Función de verosimilitud Bernoulli
Retomemos el ejemplo del experimento Bernoulli, iniciamos con una función de distribución que describa nuestro conocimiento inicial \(p(\theta)\) tal que podamos calcular la distribución objetivo \(P(\theta) = \mathcal{L}(\theta)p(\theta)\) con facilidad. En este caso elegimos una densidad beta y podemos usar función dbeta de R:
# p(theta) con theta = 0.4, a = 2, b = 2
dbeta(0.4, 2, 2)
#> [1] 1.44
# Definimos la distribución inicial
prior <- function(a = 1, b = 1){
function(theta) dbeta(theta, a, b)
}
También necesitamos especificar la función de verosimilitud, en nuestro caso tenemos repeticiones de un experimento Bernoulli por lo que: \[\mathcal{L}(\theta) = \theta^{z}(1-\theta)^{N-z}\]
y en R:
# Verosimilitid binomial
likeBern <- function(z, N){
function(theta){
theta ^ z * (1 - theta) ^ (N - z)
}
}
La distribución posterior \(p(\theta|x)\) es, por la regla de Bayes, proporcional a \(p(x|\theta)p(\theta)\). Usamos este producto como la distribución objetivo \(P(\theta)\) en el algoritmo de Metrópolis.
# posterior no normalizada
postRelProb <- function(theta){
mi_like(theta) * mi_prior(theta)
}
Implementemos el algoritmo con una inicial Beta(1,1) (uniforme) y observaciones \(z = \sum{x_i}=11\) y \(N = 14\), es decir lanzamos 14 volados de los cuales 11 resultan en éxito (águila).
# Datos observados
N = 14
z = 11
# Defino mi inicial y la verosimilitud
mi_prior <- prior() # inicial uniforme
mi_like <- likeBern(z, N) # verosimilitud de los datos observados
Proponemos una distribución normal \(\mathcal{N}(0,0.1)\) como distribución propuesta:
# para cada paso decidimos el movimiento de acuerdo a la siguiente función
caminaAleat <- function(theta){ # theta: valor actual
salto_prop <- rnorm(1, 0, sd = 0.1) # salto propuesto
theta_prop <- theta + salto_prop # theta propuesta
if(theta_prop < 0 | theta_prop > 1){ # si el salto implica salir del dominio
return(theta)
}
u <- runif(1)
p_move <- min(postRelProb(theta_prop) / postRelProb(theta), 1) # prob mover
if(p_move > u){
return(theta_prop) # aceptar valor propuesto
}
else{
return(theta) # rechazar
}
}
# Realizamos 6,000 pasos
set.seed(47405)
pasos <- 6000
camino <- numeric(pasos) # vector que guardará las simulaciones
camino[1] <- 0.1 # valor inicial
# Generamos la caminata aleatoria
for (j in 2:pasos){
camino[j] <- caminaAleat(camino[j - 1])
}
caminata <- data.frame(pasos = 1:pasos, theta = camino)
# Graficamos los primeros 3,000 pasos
ggplot(caminata[1:3000, ], aes(x = pasos, y = theta)) +
geom_point(size = 0.8) +
geom_path(alpha = 0.5) +
scale_y_continuous(expression(theta), limits = c(0, 1)) +
scale_x_continuous("Tiempo") +
geom_vline(xintercept = 600, color = "red", alpha = 0.5)
# Excluímos las primeras observaciones (etapa de calentamiento)
caminata_f <- filter(caminata, pasos > 600)
ggplot(caminata_f, aes(x = theta)) +
geom_density(adjust = 2, aes(color = "posterior")) +
labs(title = "Distribución posterior",
y = expression(p(theta)),
x = expression(theta)) +
stat_function(fun = mi_prior, aes(color = "inicial")) + # inicial
xlim(0, 1)
Si la distribución objetivo es muy dispersa y la distribución propuesta muy estrecha, entonces se necesitarán muchos pasos para que la caminata aleatoria cubra la distribución objetivo con una muestra representativa.
Por otra parte, si la distribución propuesta es muy dispersa podemos caer en rechazar demasiados valores propuestos. Imaginemos que \(\theta_{actual}\) se ubica en una zona de densidad alta, entonces cuando los valores propuestos están lejos del valor actual se ubicarán en zonas de menor densidad y \(p(\theta_{propuesta})/p(\theta_{actual})\) tenderá a ser chico y el movimiento propuesto será aceptado en pocas ocasiones.
Ejercicio:
¿Qué porcentaje de los valores propuestos fueron aceptados en el ejemplo anterior?
Si cambia la desviación estándar de la distribución propuesta a \(\sigma = 0.01\) y \(\sigma = 2\), ¿Cuál cambia el porcentaje de aceptación de valores propuestos en cada uno de los dos casos?
¿De los 3 valores que usamos, qué desviación estándar crees que sea más conveniente?
mean(caminata_f$theta)
#> [1] 0.748
sd(caminata_f$theta)
#> [1] 0.109
sims_y <- rbinom(nrow(caminata_f), size = 1, prob = caminata_f$theta)
mean(sims_y) # p(y = 1 | x) probabilidad predictiva
#> [1] 0.75
sd(sims_y)
#> [1] 0.433
Consideramos la situación en la que nos interesa estudiar dos proporciones \(\theta_1\) y \(\theta_2\) correspondientes a dos grupos. Queremos determinar que debemos creer de estas proporciones tras observar datos provenientes de ambos grupos.
En el marco Bayesiano comenzamos definiendo nuestras creencias iniciales. En este caso nuestras creencias describen combinaciones de parámetros, es decir debemos especificar la distribución inicial bivariada \(p(\theta_1, \theta_2)\) para todas las combinaciones \(\theta_1, \theta_2\).
Un caso sencillo es asumir que nuestras creencias de \(\theta_1\) son independientes de nuestras creencias de \(\theta_2\). Esto implica:
\[p(\theta_1, \theta_2) = p(\theta_1)p(\theta_2)\] para todo valor de \(\theta_1\) y de \(\theta_2\) y donde \(p(\theta_1)\) y \(p(\theta_2)\) corresponden a las distribuciones iniciales marginales. Las creencias de los dos parámetros no tienen porque ser independientes; por ejemplo, puedo creer que dos monedas acuñadas en la misma casa de moneda tienen un sesgo similar, en este caso las manipulaciones matemáticas son más complicadas pero es posible describir estas creencias mediante una distribución conjunta bivariada.
Además de las creencias iniciales tenemos datos observados. Suponemos que los lanzamientos dentro de cada grupo son independientes y que los lanzamientos entre los grupos también lo son. Es importante recalcar que siempre suponemos independencia en los datos, sin importar nuestros supuestos de independencia en las creencias.
Para el primer grupo observamos la secuencia \(\{x_{1,1},...,x_{1,N_1}\}\) que contiene \(z_1\) águilas, y en el otro grupo observamos la sucesión \(\{x_{2,1},...,x_{2,N_2}\}\) que contiene \(z_2\) águilas. Es decir, \[z_1 = \sum_{i=1}^{N_1}x_{1,i}, \quad z_2 = \sum_{i=1}^{N_2}x_{2,i}\] Por simplicidad denotamos todos los datos por \(x =\{x_{1,1},...,x_{1,N_1}x_{2,1},...,x_{2,N_2}\}\)
Debido a la independencia de los lanzamientos tenemos la siguiente función de verosimilitud:
\[p(x|\theta_1,\theta_2)=\prod_{i=1}^{N_1}p(x_{1,i}|\theta_1,\theta_2)\cdot \prod_{i=1}^{N_2}p(x_{2,i}|\theta_1,\theta_2)\] \[= \theta_1^{z_1}(1-\theta_1)^{N_1-z_1}\cdot \theta_2^{z_2}(1-\theta_2)^{N_2-z_2}\]
Usamos la regla de Bayes para calcular la distribución posterior:
\[p(\theta_1,\theta_2|x)=p(x|\theta_1,\theta_2)p(\theta_1, \theta_2) / p(x)\] \[= \frac{p(x|\theta_1,\theta_2)p(\theta_1, \theta_2)} { \int\int p(x|\theta_1,\theta_2)p(\theta_1, \theta_2)d\theta_1d\theta_2}\] Resta calcular la doble integral en el denominador que define la evidencia \(p(x)\).
Siguiendo el caso de la familia Beta-Bernoulli univariada que estudiamos en el caso de una proporción y suponiendo independencia en las creencias, \(p(\theta_1, \theta_2) = p(\theta_1)p(\theta_2)\). Escribimos la densidad inicial como el producto de dos densidades beta, donde \(\theta_1\) se distribuye \(Beta(a_1, b_1)\) y \(\theta_2\) se distribuye \(Beta(a_2, b_2)\).
\[p(\theta_1, \theta_2)=\frac{\theta_1^{a_1-1}(1-\theta)^{b_1-1}} {B(a_1, b_1)}\frac{ \theta_2^{a_2-1}(1-\theta)^{b_2-1}}{B(a_2, b_2)}\]
donde la función beta está definida en términos de la función gamma de Euler \(B(a, b) = \Gamma(a)\Gamma(b) / \Gamma(a+b)\).
La posterior se escribe:
\[p(\theta_1, \theta_2|x)=\frac{\theta_1^{z_1+a_1-1}(1-\theta)^{N_1-z_1+b_1-1}\cdot \theta_2^{z_2+a_2-1}(1-\theta)^{N_2-z_2+b_2-1}}{p(x)\cdot B(a_1, b_1) \cdot B(a_2, b_2)}\] Vemos que, cuando la inicial es el producto de distribuciones beta independientes, la posterior también es el producto de distribuciones beta independientes. Por lo tanto
\[p(x) = \frac{B(z_1+a_1,N_1-z_1+b_1)B(z_2,N_2-z_2+_b2)}{B(a_1, b_1) \cdot B(a_2, b_2)}\] Visualizaremos las distribuciónes bivariadas mediante gráficas de contornos.
Consideremos la siguinete elección de los parámetros iniciales \(a_1=b_1=a_2=b_2=3\). Por lo que la distribución inicial bivariada es:
grid <- expand.grid(x = seq(0.01, 1, 0.01), y = seq(0.01, 1, 0.01))
grid_inicial <- grid %>%
mutate(z = dbeta(x, 3, 3) * dbeta(y, 3, 3))
binom_2_inicial <- ggplot(grid_inicial, aes(x = x, y = y, z = z)) +
# geom_tile(aes(fill = z)) +
geom_raster(aes(fill = z), show.legend = FALSE) +
geom_contour(colour = "white") +
# stat_contour(binwidth = 0.6, aes(color = ..level..), show.legend = FALSE) +
scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
scale_y_continuous(expression(theta[2]), limits = c(0, 1)) +
scale_color_gradient(expression(p(theta[1],theta[2])), limits = c(0, 8.6)) +
scale_fill_gradient(expression(Probabilidad), limits = c(0, 8.6)) +
coord_fixed() +
labs(title = "Inicial")
Si asumimos que los datos se pueden resumir mediante \(z_1=5, N_1=7, z_2=2, N_2=7\), la función de verosimiliud es:
grid_v <- grid %>%
mutate(z = dbeta(x, 6, 3) * dbeta(y, 3, 6))
binom_2_verosimilitud <- ggplot(grid_v, aes(x = x, y = y, z = z)) +
geom_raster(aes(fill = z), show.legend = FALSE) +
geom_contour(colour = "white") +
# stat_contour(binwidth = 0.6, aes(color = ..level..), show.legend = FALSE) +
scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
scale_y_continuous(expression(theta[2]), limits = c(0, 1)) +
scale_color_gradient(expression(p(theta[1],theta[2])), limits = c(0, 8.6)) +
scale_fill_gradient(expression(Probabilidad), limits = c(0, 8.6)) +
coord_fixed() +
labs(title = "Verosimilitud")
Note que usamos un producto de distribuciones beta para expresar el producto de distribuciones Bernoulli (la normalización no coincide).
La distribución posterior bivariada es entonces:
grid_post <- grid %>%
mutate(z = dbeta(x, 8, 5) * dbeta(y, 5, 8))
binom_2_posterior <- ggplot(grid_post, aes(x = x, y = y, z = z)) +
geom_raster(aes(fill = z)) +
geom_contour(colour = "white", show.legend = FALSE) +
# stat_contour(binwidth = 0.6, aes(color = ..level..), show.legend = FALSE) +
scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
scale_y_continuous(expression(theta[2]), limits = c(0, 1)) +
# scale_color_gradient(expression(p(theta[1],theta[2])), limits = c(0, 8.6)) +
scale_fill_gradient(expression(Probabilidad), limits = c(0, 8.6)) +
coord_fixed() +
labs(title = "Posterior")
grid.arrange(binom_2_inicial, binom_2_verosimilitud, binom_2_posterior,
nrow = 1, widths = c(0.3, 0.3, 0.4))
También podemos realizar gráficas tridimensionales:
# library(plotly)
# grid_inicial_pl <- grid_inicial %>% spread(y, z) %>% as.matrix()
# pl_inicial <- plot_ly(z = grid_inicial_pl) %>% add_surface(cmin = 0, cmax = 9)
# grid_v_pl <- grid_v %>% spread(y, z) %>% as.matrix()
# pl_verosimilitud <- plot_ly(z = grid_v_pl) %>% add_surface(cmin = 0, cmax = 9)
# grid_post_pl <- grid_post %>% spread(y, z) %>% as.matrix()
# pl_post <- plot_ly(z = grid_post_pl) %>% add_surface(cmin = 0, cmax = 9)
#
# pl_inicial
# pl_verosimilitud
# pl_post
Al igual que en el caso univariado escribimos la la distribución inicial, la función de verosimilitud y la distribución posterior no normalizada:
# Inicial bivariada
prior2 <- function(a_1 = 3, b_1 = 3, a_2 = 3, b_2 = 3){
function(theta) dbeta(theta[1], a_1 , b_1) * dbeta(theta[2], a_2 , b_2)
}
# Verosimilitud binomial bivariada
likeBern2 <- function(z_1, N_1, z_2, N_2){
function(theta){
theta[1] ^ z_1 * (1 - theta[1]) ^ (N_1 - z_1) *
theta[2] ^ z_2 * (1 - theta[2]) ^ (N_2 - z_2)
}
}
# Posterior bivariada no normalizada
postRelProb2 <- function(theta){
mi_like(theta) * mi_prior(theta)
}
Implementemos el algoritmo con el producto de dos distribuciones iniciales Beta(3,3) y observaciones \(z_1=5\), \(N = 7\), \(z_1=2\) y \(N = 7\), es decir lanzamos 14 volados, 7 de cada moneda.
z_1=5; N_1=7; z_2=2; N_2=7; a_1=3; a_2=3; b_1=3; b_2=3
# Datos observados
N = 14
z = 11
# Defino mi inicial y la verosimilitud
mi_prior <- prior2()
mi_like <- likeBern2(z_1 = z_1, N_1 = N_1, z_2 = z_2, N_2 = N_2) # verosimilitud
Para proponer saltos usaremos una distribución propuesta normal bivariada.
Recordemos el algoritmo: el movimiento se acepta de manera definitiva si la distribución objetivo (posterior no normalizada) es más densa en la posición propuesta que en la posición actual, y se acepta de manera probabilistica si la distribución objetivo es más densa en posición actual que en la posición propuesta.
# para cada paso decidimos el movimiento de acuerdo a la siguiente función
caminaAleat2 <- function(theta){ # theta: valor actual
salto_prop <- MASS::mvrnorm(n = 1 , mu = rep(0, 2),
Sigma = matrix(c(0.1, 0, 0, 0.1), byrow = TRUE, nrow = 2)) # salto propuesto
theta_prop <- theta + salto_prop # theta propuesta
if(all(theta_prop < 0) | all(theta_prop > 1)){ # salir del dominio
return(theta)
}
u <- runif(1)
p_move = min(postRelProb2(theta_prop) / postRelProb2(theta), 1) # prob mover
if(p_move > u){
return(theta_prop) # aceptar valor propuesto
}
else{
return(theta) # rechazar
}
}
set.seed(47405)
pasos <- 6000
camino <- matrix(0, nrow = pasos, ncol = 2) # vector que guardará las simulaciones
camino[1, ] <- c(0.1, 0.8) # valor inicial
# Generamos la caminata aleatoria
for (j in 2:pasos){
camino[j, ] <- caminaAleat2(camino[j - 1, ])
}
caminata <- data.frame(pasos = 1:pasos, theta_1 = camino[, 1],
theta_2 = camino[, 2])
Graficamos la caminata aleatoria:
ggplot(caminata, aes(x = theta_1, y = theta_2)) +
geom_point(size = 0.8) +
geom_path(alpha = 0.3) +
scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
scale_y_continuous(expression(theta[2]), limits = c(0, 1)) +
coord_fixed()
También podemos graficar las caminatas aleatorias “marginales”:
caminata_m <- caminata %>%
gather(parametro, val, theta_1, theta_2) %>%
arrange(pasos)
ggplot(caminata_m[1:2000, ], aes (x = pasos, y = val)) +
geom_path(alpha = 0.5) +
facet_wrap(~parametro, ncol = 1) +
scale_y_continuous("", limits = c(0, 1))
El algoritmo de Metrópolis es muy general y se puede aplicar a una gran variedad de problemas. Sin embargo, afinar los parámetros de la distribución propuesta para que el algoritmo funcione correctamente, es decir, para que explore de manera eficiente el espacio de parámetros de la distribucion objetivo puede ser complicado. Por otra parte, como veremos a continuación, el muestredor de Gibbs no necesita de una distribución propuesta, ni de la distribución objetivo conjunta.
Para implementar un muestreador de Gibbs se necesita únicamente ser capaz de generar muestras de las distribuciones posteriores marginales de cada parámetro individual condicinales al resto de los parámetros. Esto es, el muestreador de Gibbs permite generar muestras de la distribucion posterior: \[p(\theta_1,...,\theta_p|x)\] siempre y cuando podamos generar valores de todas las distribuciones posteriores marginales condicionales: \[p(\theta_k|\theta_1,...,\theta_{k-1},\theta_{k+1},...,\theta_p,x)\]
El proceso del muestreador de Gibbs es una caminata aleatoria a lo largo del espacio de parámetros. La caminata inicia en un punto arbitrario y, en cada tiempo, el siguiente paso depende únicamente de la posición actual. Por tanto, el muestredor de Gibbs es un proceso Cadena de Markov vía Monte Carlo (MCMC). La diferencia entre el muestreador de Gibbs y el algoritmo de Metrópolis radica en cómo se eligen los siguientes pasos y en el hecho de que un nuevo paso siempre se acepta.
En el caso del muestreador de Gibbs, en cada punto de la caminata aleatoria se selecciona uno de los componentes del vector de parámetros (típicamente se cicla en orden). Supongamos que se selecciona el parámetro \(\theta_k\), entonces obtenemos un nuevo valor para este parámetro generando una simulación de la distribución condicional \[p(\theta_k|\theta_1,...,\theta_{k-1},\theta_{k+1},...,\theta_p,x)\]
El nuevo valor \(\theta_k\) junto con los valores que aún no cambian \(\theta_1,...,\theta_{k-1},\theta_{k+1},...,\theta_p\) constituyen la nueva posición en la caminata aleatoria.
Seleccionamos una nueva componente del vector de parámetros \(\theta\) y repetimos el proceso.
El muestreador de Gibbs es útil cuando no podemos determinar de manera analítica la distribución posterior conjunta (no normalizada) pero sí podemos determinar todas las distribuciones posteriores marginales condicionales y simular de ellas.
Ejemplo: Ejemplificaremos el muestreador de Gibbs con el ejemplo de las dos proporciones binomiales, a pesar de no ser necesario en este caso.
Comenzamos identificando las distribuciones condicionales posteriores para cada parámetro:
\[p(\theta_1|\theta_2,x) = p(\theta_1,\theta_2|x) / p(\theta_2|x)\] \[= \frac{p(\theta_1,\theta_2|x)} {\int p(\theta_1,\theta_2|x) d\theta_1}\]
Usando iniciales \(Beta(a_1, b_1)\) y \(Beta(a_2,b_2)\), obtenemos:
\[p(\theta_1|\theta_2,x) = \frac{Beta(\theta_1|z_1 + a_1, N_1 - z_1 + b_1) Beta(\theta_2|z_2 + a_2, N_2 - z_2 + b_2)}{\int Beta(\theta_1|z_1 + a_1, N_1 - z_1 + b_1) Beta(\theta_2|z_2 + a_2, N_2 - z_2 + b_2) d\theta_1}\] \[= \frac{Beta(\theta_1|z_1 + a_1, N_1 - z_1 + b_1) Beta(\theta_2|z_2 + a_2, N_2 - z_2 + b_2)}{Beta(\theta_2|z_2 + a_2, N_2 - z_2 + b_2)}\] \[=Beta(\theta_1|z_1 + a_1, N_1 - z_1 + b_1)\]
Debido a que la posterior es el producto de dos distribuciones Beta independientes es claro que \(p(\theta_1|\theta_2,x)=p(\theta_1|x)\).
Una vez que determinamos las distribuciones condicionales, simplemente hay que encontrar una manera de obtener muestras de éstas, en R podemos usar la función rbeta.
pasos <- 12000
camino <- matrix(0, nrow = pasos, ncol = 2) # vector que guardará las simulaciones
camino[1, 1] <- 0.1 # valor inicial
camino[1, 2] <- 0.1
# Generamos la caminata aleatoria
for (j in 2:pasos){
if(j %% 2){
camino[j, 1] <- rbeta(1, z_1 + a_1, N_1 - z_1 + b_1)
camino[j, 2] <- camino[j - 1, 2]
}
else{
camino[j, 2] <- rbeta(1, z_2 + a_2, N_2 - z_2 + b_2)
camino[j, 1] <- camino[j - 1, 1]
}
}
caminata <- data.frame(pasos = 1:pasos, theta_1 = camino[, 1],
theta_2 = camino[, 2])
Graficamos la caminata aleatoria bivariada:
ggplot(caminata[1:2000, ], aes(x = theta_1, y = theta_2)) +
geom_point(size = 0.8) +
geom_path(alpha = 0.5) +
scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
scale_y_continuous(expression(theta[2]), limits = c(0, 1)) +
coord_fixed()
También podemos graficar las caminatas aleatorias marginales:
caminata_g <- filter(caminata, pasos %% 2 == 0) %>%
gather(parametro, val, theta_1, theta_2) %>%
mutate(pasos = rep(1:6000, 2)) %>%
arrange(pasos)
ggplot(caminata_g[1:2000, ], aes(x = pasos, y = val)) +
geom_path(alpha = 0.3) +
facet_wrap(~parametro, ncol = 1) +
scale_y_continuous("", limits = c(0, 1))
# Metropolis
caminata_m %>%
filter(pasos > 1000) %>% # eliminamos el calentamiento
group_by(parametro) %>%
summarise(
media = mean(val),
mediana = median(val),
std = sd(val)
)
#> # A tibble: 2 x 4
#> parametro media mediana std
#> <chr> <dbl> <dbl> <dbl>
#> 1 theta_1 0.605 0.606 0.124
#> 2 theta_2 0.392 0.388 0.136
# Gibbs
caminata_g %>%
filter(pasos > 1000) %>%
group_by(parametro) %>%
summarise(
media = mean(val),
mediana = median(val),
std = sd(val)
)
#> # A tibble: 2 x 4
#> parametro media mediana std
#> <chr> <dbl> <dbl> <dbl>
#> 1 theta_1 0.619 0.630 0.129
#> 2 theta_2 0.384 0.378 0.131
Note que en ambos casos las estimaciones son muy cercanas.
También podemos comparar los sesgos de las dos monedas:
caminata <- caminata %>%
mutate(dif = theta_1 - theta_2)
ggplot(caminata, aes(x = dif)) +
geom_histogram(fill = "gray") +
geom_vline(xintercept = 0, alpha = 0.8, color = "red")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
La principal ventaja del muestreador de Gibbs sobre el algoritmo de Metrópolis es que no hay necesidad de seleccionar una distribución propuesta y no hay que lidiar con lo ineficiente de rechazar valores, además de que no necesitamos conocer la distrbución objetivo conjunta. A cambio, debemos ser capaces de derivar las distribuciones posteriores marginales condicionales de cada parámetro y generar muestras de estas.
Retomemos el caso de observaciones normales. Supongamos que tenemos una muestra \(x_1,...,x_N\) de observaciones independientes e identicamente distribuidas i.i.d., con \(x_i \sim N(\mu, \sigma^2)\). Vamos a considerar el caso de media desconocida, el caso de varianza desconocida y el caso en el que ambas son desconocidas.
Normal con media desconocida. Supongamos que \(\sigma^2\) es conocida, por lo que nuestro parámetro de interés es únicamente \(\mu\). Describiremos nuestro conocimiento inicial de \(\mu\) a través de una distribución inicial normal con media \(m\) y varianza \(\tau^2\) conocidas, ya que ésta es una distribución conjugada de la verosimilitud en cuestión: \[\mu \sim N(m, \tau^2)\] La distribución posterior resulta en: \[\mu|x \sim N\bigg(\frac{\sigma^2}{\sigma^2 + N\tau^2}m + \frac{N\tau^2}{\sigma^2 + N \tau^2}\bar{x}, \frac{\sigma^2 \tau^2}{\sigma^2 + N\tau^2}\bigg)\]
Normal con varianza desconocida. Supongamos que \(\mu\) es conocida, por lo que nuestro parámetro de interés es únicamente \(\sigma^2\). En este caso una distribución conveniente para describir nuestro conocimiento inicial es la distribución Gamma Inversa con parámetros conocidos, ya que ésta es conjugada a la verosimilitud en cuestión.
La distribución Gamma Inversa es una distribución contínua con dos parámetros y que toma valores en los positivos. Como su nombre lo indica, esta distribución corresponde al recírpoco de una variable cuya distribución es Gamma. Recordemos que si \(x\sim Gamma(\alpha, \beta)\) entonces:
\[p(x)=\frac{1}{\beta^{\alpha}\Gamma(\alpha)}x^{\alpha-1}e^{-x/\beta}\]
donde \(x>0\). Ahora si \(y\) es la variable aleatoria recírpoca de \(x\) entonces:
\[p(y)=\frac{\beta^\alpha}{\Gamma(\alpha)}y^{-\alpha - 1} exp{-\beta/y}\]
con media \[\frac{\beta}{\alpha-1}\] y varianza \[\frac{\beta^2}{(\alpha-1)^2(\alpha-2)}.\]
Debido a la relación entre las distribuciones Gamma y Gamma Inversa, podemos utilizar la función rgamma de R para generar valores de distribución gamma inversa.
# 1. simulamos valores porvenientes de una distribución gamma
x_gamma <- rgamma(2000, shape = 5, rate = 1)
# 2. invertimos los valores simulados
x_igamma <- 1 / x_gamma
# También podemos usar las funciones de MCMCpack
library(MCMCpack)
#> Loading required package: coda
#> Loading required package: MASS
#>
#> Attaching package: 'MASS'
#> The following object is masked from 'package:plotly':
#>
#> select
#> The following object is masked from 'package:dplyr':
#>
#> select
#> ##
#> ## Markov Chain Monte Carlo Package (MCMCpack)
#> ## Copyright (C) 2003-2018 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
#> ##
#> ## Support provided by the U.S. National Science Foundation
#> ## (Grants SES-0350646 and SES-0350613)
#> ##
#>
#> Attaching package: 'MCMCpack'
#> The following object is masked from 'package:LearnBayes':
#>
#> rdirichlet
x_igamma <- data.frame(x_igamma)
ggplot(x_igamma, aes(x = x_igamma)) +
geom_histogram(aes(y = ..density..), binwidth = 0.05, fill = "gray") +
stat_function(fun = dinvgamma, args = list(shape = 5, scale = 1),
color = "red")
Volviendo al problema de inferir acerca del parámetro \(\sigma^2\), si resumimos nuestro conocimiento inicial a través de una distribución Gamma Inversa tenemos \[p(\sigma^2)=\frac{\beta^\alpha}{\Gamma(\alpha)}\frac{1}{(\sigma^2)^{\alpha + 1}} e^{-\beta/\sigma^2}\]
Recordemos la función de verosimiltud: \[p(x|\mu, \sigma^2)=\frac{1}{(2\pi\sigma^2)^{N/2}}exp\left(-\frac{1}{2\sigma^2}\sum_{j=1}^{N}(x_j-\mu)^2\right)\]
La distribución posterior resulta en:
\[p(\sigma^2) \propto p(x|\mu,\sigma^2)p(\sigma^2)\]
por lo que obtenemos una distribución Gamma Inversa \[\sigma^2|x \sim GI(N/2+\alpha, \beta + 1/2 \sum(x_i - \mu)^2)\].
Por tanto, como anticipamos, tenemos que la inicial Gamma Inversa con verosimilitud Normal para el parámetro \(\sigma^2\) es una familia conjugada.
Normal con media y varianza desconocidas. Sea \(\theta=(\mu, \sigma^2)\). Especificamos la siguiente inicial para \(\theta\): \[\theta \sim N(\mu|m, \tau^2)\cdot GI(\sigma^2|\alpha, \beta)\] donde suponemos que los hiperparámetros \(m,\tau^2, \alpha, \beta\) son conocidos.
Entonces, la distribución posterior no normalizada es: \[ p(\theta|x) \propto p(x|\theta) p(\theta)\] \[= \frac{1}{(\sigma^2)^{N/2}} exp\bigg(-\frac{1}{2\sigma^2}\sum_{i=1}^N (x_i-\mu)^2 \bigg) exp\bigg(-\frac{1}{2\tau^2}(\mu-m)^2)\bigg) \frac{1}{(\sigma^2)^{\alpha +1}} exp\bigg(-\frac{\beta}{\sigma^2}\bigg)\]
En esta última distribución no reconocemos el núcleo de niguna distribución conocida, por lo que no podemos utilizar el método de distribuciones conjugadas.
Sin embargo, si nos concentramos únicamente en los términos que involucran a \(\mu\) tenemos:
\[exp\left(-\frac{1}{2}\left( \mu^2 \left( \frac{N}{\sigma^2} + \frac{1}{\tau^2} \right) - 2\mu\left(\frac{\sum_{i= 1}^n x_i}{\sigma^2} + \frac{m}{\tau^2}\right) \right)\right)\]
Esta expresión depende de \(\mu\) y \(\sigma^2\), por lo que condicional a \(\sigma^2\) observamos el núcleo de una distribución normal:
\[\mu|\sigma^2,x \sim N\left(\frac{n\tau^2}{n\tau^2 + \sigma^2}\bar{x} + \frac{\sigma^2}{N\tau^2 + \sigma^2}m, \frac{\tau^2\sigma^2}{n\tau^2 + \sigma^2} \right)\]
Si nos fijamos únicamente en los tárminos que involucran a \(\sigma^2\) tenemos:
\[\frac{1}{(\sigma^2)^{N/2+\alpha+1}}exp\left(- \frac{1}{\sigma^2} \left(\sum_{i=1}^N \frac{(x_i-\mu)^2}{2} + \beta \right) \right)\]
por lo que condicional a \(\mu\) reconocemos el núcleo de una distribución Gamma Inversa:
\[\sigma^2|\mu,x \sim GI\left(\frac{N}{2} + \alpha, \sum_{i=1}^n \frac{(x_i-\mu)^2}{2} + \beta \right)\]
Obtenemos de esta forma las distribuciones posteriores marginales condicionales completas \(p(\mu|\sigma^2, x)\) y \(p(\sigma^2|\mu, x)\), las cuales son distribuciones conocidas y de las cuales podemos simular.
Ésto nos permite implementar un muestreador de Gibbs.
Comenzamos definiendo las distrbuciones iniciales:
\(\mu \sim N(1.5, 16)\), esto es \(m = 1.5\) y \(\tau^2 = 16\).
\(\sigma^2 \sim GI(3, 3)\), esto es \(\alpha = \beta = 3\).
Ahora supongamos que observamos 50 realizaciones (observaciones) provenientes de la distribución de interés (la que da lugar a la función de verosimilitud):
N <- 50 # Observamos 50 realizaciones
set.seed(122)
x <- rnorm(N, 2, 2)
x
#> [1] 4.6214 0.2483 2.3990 2.9319 -1.6041 4.8975 2.5977 2.7236
#> [9] -0.0139 1.4860 1.7357 0.3167 2.5485 -2.9252 -2.3068 4.3184
#> [17] 3.3795 3.7605 0.1133 3.4381 0.9243 0.9547 -0.1058 2.2030
#> [25] 5.7270 1.9608 -0.1566 2.3452 3.0661 5.9045 4.8227 3.2027
#> [33] 0.1720 5.1609 1.0602 5.2037 2.7455 2.0678 2.2082 -2.0367
#> [41] 0.6840 -1.2170 -0.6882 1.6067 4.7513 2.3466 1.1214 -1.1906
#> [49] 0.5114 5.6304
Escribimos el muestreador de Gibbs:
m <- 1.5; tau2 <- 16; alpha <- 3; beta <- 3 # parámetros de iniciales
pasos <- 20000
camino <- matrix(0, nrow = pasos + 1, ncol = 2) # vector que guardará las simulaciones
camino[1, 1] <- 0 # valor inicial media
camino[1, 2] <- 0.1 # valor inicial varianza
# Generamos la caminata aleatoria
for (j in 2:(pasos)){
# sigma^2
if(j %% 2){
mu <- camino[j - 1, 1] #Actualizar los par?metros
a <- N / 2 + alpha
b <- sum((x - mu) ^ 2) / 2 + beta
camino[j, 1] <- camino[j - 1, 1] #mu no cambia
camino[j, 2] <- 1/rgamma(1, shape = a, rate = b) #Actualizar sigma2
}
# mu
else{
sigma2 <- camino[j - 1, 2] #Actualizar los par?metros
media <- (N * tau2 * mean(x) + sigma2 * m) / (N * tau2 + sigma2)
var <- sigma2 * tau2 / (N * tau2 + sigma2)
camino[j, 2] <- camino[j - 1, 2] #sigma2 no cambia
camino[j, 1] <- rnorm(1, media, sd = sqrt(var)) #Actualizar mu
}
}
caminata <- data.frame(pasos = 1:pasos, mu = camino[1:pasos, 1],
sigma2 = camino[1:pasos, 2])
Podemos graficar la caminata aleatoria:
ggplot(caminata[1000:2000, ], aes(x = mu, y = sigma2)) +
geom_point(size = 0.8) +
geom_path(alpha = 0.5) +
scale_x_continuous(expression(mu), limits = c(0, 5)) +
scale_y_continuous(expression(sigma^2), limits = c(0, 10)) +
coord_fixed()
#> Warning: Removed 4 rows containing missing values (geom_point).
También podemos graficar las caminatas aleatorias marginales (eliminando el periodo de calentamiento):
caminata_g <- caminata %>%
gather(parametro, val, mu, sigma2) %>%
arrange(pasos)
ggplot(filter(caminata_g, pasos > 15000), aes(x = pasos, y = val)) +
geom_path(alpha = 0.3) +
facet_wrap(~parametro, ncol = 1, scales = "free") +
scale_y_continuous("")
ggplot(filter(caminata_g, pasos > 5000), aes(x = val)) +
geom_histogram(fill = "gray") +
facet_wrap(~parametro, ncol = 1, scales = "free")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
y calcular algunos resúmenes de la posterior:
caminata_g %>%
filter(pasos > 1000) %>% # eliminamos la etapa de calentamiento
group_by(parametro) %>%
summarise(
mean(val),
sd(val),
median(val)
)
#> # A tibble: 2 x 4
#> parametro `mean(val)` `sd(val)` `median(val)`
#> <chr> <dbl> <dbl> <dbl>
#> 1 mu 1.91 0.305 1.91
#> 2 sigma2 4.75 0.953 4.61
\[p(y) =\int p(y|\theta)p(\theta|x)d\theta\]
Por tanto podemos aproximar la distribución predictiva posterior como:
caminata_f <- filter(caminata, pasos > 5000)
caminata_f$y_sims <- rnorm(1:nrow(caminata_f), caminata_f$mu, caminata_f$sigma2)
ggplot(caminata_f, aes(x = y_sims)) +
geom_histogram(fill = "gray") +
geom_vline(aes(xintercept = mean(y_sims)), color = "red")
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Ejercicio: ¿Cuál es la probabilidad de que una observación futura sea mayor a 8?
Ejercicio: En estadística bayesiana es común parametrizar la distribución Normal en términos de la precisión \(\nu\), la cual se define como el inverso de la varianza \(\sigma^2\). Si parametrizamos de esta manera \(\nu = 1/\sigma^2\), repita el ejemplo anterior con la diferencia de utilizar la distribución inicial Gamma para \(\nu\) en lugar de la Gamma inversa para \(\sigma^2\).
JAGS (Just Another Gibbs Sampler), WinBUGS y OpenBUGS son programas que implementan métodos MCMC para generar simulaciones de distribuciones posteriores. Los paquetes rjags y R2jags permiten ajustar modelos en JAGS desde R. Es muy fácil utilizar estos programas pues uno simplemente debe especificar las distribuciones iniciales, la verosimilitud y los datos observados.
Instalar JAGS.
Instalar los paquetes R2jags y rjags de R.
Ilustraremos el proceso de simulación en JAGS para el caso del sesgo de una moneda \(\theta\) que consideramos anteriormente.
Primero, especificamos los datos observados:
# N = 14; z = 11
N <- 14
x <- c(rep(0, 3), rep(1, 11)) #Datos simulados
z <- sum(x)
z
#> [1] 11
Consideremos el modelo gráfico asociado:
library(DiagrammeR)
grViz("
digraph betabinom {
# a 'graph' statement
graph [overlap = true]
node [shape = box,penwidth = 0.1, fixedsize = true, fontsize = 9,
fontname = Helvetica, width = 0.2, height = 0.2]
ab[label = 'a,b'];
node [shape = circle,
fixedsize = true, fontsize = 9,
fontname = Helvetica, width = 0.2] // sets as circles
theta[label = 'θ']; y;
# several 'edge' statements
edge[color = grey, arrowsize = 0.5, penwidth = 0.5]
ab->theta; theta->y;
}
")
donde \(a,b\) son los hiperparámetros de la distribución inicial Beta de \(\theta\), los cuales los elegimos de tal forma que representen nuestras creencias iniciales acerca de \(\theta\).
El diagrama captura las dependencias entre los datos \(y\), el parámetro \(\theta\) y los hiperparámetros \(a,b\). Ésto facilita la implementación en JAGS pues cada flecha en el diagrama corresponde a una línea de código en la especificación del modelo.
Especifiquemos el modelo:
modelo_bb.bugs <-
'
model{
for(i in 1:N){
x[i] ~ dbern(theta)
}
# inicial
theta ~ dbeta(1, 1)
}
'
El ciclo for indica que cada dato observado \(x_i\) proviene de una distribución Bernoulli con parámetro \(\theta\). Afuera del ciclo escribimos la distribución inicial, \(\theta \sim Beta(1, 1)\).
El modelo ya esta especificado. Aún hace falta especificar el valor inicial de los parámetros \(\theta\). JAGS tiene una manera de hacerlo automaticamente, pero muchas veces vale la pena tener control de los valores iniciales. En ocasiones la eficiencia del proceso puede incrementar si seleccionamos valores iniciales razonables. Kruschke sugiere utilizar como valores iniciales los estimadores de máxima verosimilitud de los parámetros; esto es porque usualmente la distribución posterior no está muy lejana de la función de verosimilitud. En este caso el estimador de máxima verosimilitud para \(\theta\) es la proporción observada \(\hat{\theta}=z/N\).
theta_init <- sum(x) / N
theta_init
#> [1] 0.786
Esta manera de especificar los valores iniciales no siempre es la más recomendada, ya que cuando queremos evaluar la convergencia de la cadena de Markov muchas veces se sugiere correr varias caminatas aleatorias con puntos iniciales dispersos a lo largo del espacio de parámetros, de tal manera que cuando las cadenas de Markov convergen se pueda determinar que la etapa de calentamiento a terminado.
Un punto medio es iniciar la caminata aleatoria en un punto aleatorio cercano al estimador de máxima verosimilitud:
init_theta <- function(){
x_s <- sample(x, replace = TRUE)
return(list(theta = sum(x_s) / N))
}
init_theta()
#> $theta
#> [1] 0.857
Ahora llamamos a JAGS y generamos las cadenas de Markov. Para esto usaremos el paquete R2jags; otro paquete para llamar JAGS desde R es rjags.
library(R2jags)
#> Warning: package 'R2jags' was built under R version 3.5.1
#> Loading required package: rjags
#> Linked to JAGS 4.3.0
#> Loaded modules: basemod,bugs
#>
#> Attaching package: 'R2jags'
#> The following object is masked from 'package:coda':
#>
#> traceplot
cat(modelo_bb.bugs, file = 'modelo_bb.bugs')
jags_fit <- jags(
model.file = "modelo_bb.bugs", # modelo de JAGS
inits = init_theta, # valores iniciales
data = list(x = x, N = N), # lista con los datos
parameters.to.save = c("theta"), # parámetros por guardar
n.chains = 1, # número de cadenas
n.iter = 1000, # número de pasos
n.burnin = 500 # calentamiento de la cadena
)
#> module glm loaded
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 14
#> Unobserved stochastic nodes: 1
#> Total graph size: 17
#>
#> Initializing model
# plot(jags_fit)
Podemos ver un resumen del ajuste:
head(jags_fit$BUGSoutput$summary)
#> mean sd 2.5% 25% 50% 75% 97.5%
#> deviance 15.522 1.406 14.55 14.659 14.96 15.811 19.036
#> theta 0.751 0.107 0.52 0.688 0.76 0.832 0.925
y graficar la cadena de Markov:
traceplot(jags_fit, varname = "theta")
Recordemos el ejemplo Normal con media y varianza desconocidas.
N <- 50 # Observamos 50 realizaciones
set.seed(122)
x <- rnorm(N, 2, 2) #Datos simulados
Ejercicio: ¿Cuál es el modelo gráfico asociado?
El modelo es:
modelo_normal.bugs <-
'
model{
for(i in 1:N){
x[i] ~ dnorm(mu, nu)
}
# iniciales
nu ~ dgamma(3, 3)
sigma2 <- 1 / nu
mu ~ dnorm(1.5, 1 / 16)
}
'
El ciclo for indica que cada dato observado \(x_i\) proviene de una distribución Normal con media \(\mu\) y precisión \(\nu\). Afuera del ciclo escribimos las distribuciones iniciales, \(\nu \sim Gamma(3, 3)\), esto es, \(\sigma^2 \sim GI(3, 3)\); y \(\mu \sim N(1.5, 16)\), esto es, \(\mu\) se distribuye Normal con media \(m=1.5\) y varianza \(\tau^2=16\). Note que en JAGS la distribución Normal se especifica en términos de la media y la precisión.
El modelo ya esta especificado, pero aun falta indicar los valores iniciales de los parámetros del modelo. Para ésto, definimos los valores en R y después los mandamos a JAGS para generar la caminata aleatoria.
library(R2jags)
cat(modelo_normal.bugs, file = 'modelo_normal.bugs')
# valores iniciales para los parámetros, si no se especifican la función jags generará valores iniciales (en este caso tomamos valores aleatorios)
jags_inits <- function(){
list(mu = rnorm(1, mean(x), 5), nu = 1 / runif(1, 2, 4))
}
jags_fit <- jags(
model.file = "modelo_normal.bugs", # modelo de JAGS
inits = jags_inits, # valores iniciales
data = list(x = x, N = N), # lista con los datos
parameters.to.save = c("mu", "sigma2"), # parámetros por guardar
n.chains = 1, # número de cadenas
n.iter = 10000, # número de pasos
n.burnin = 1000, # calentamiento de la cadena
n.thin = 1
)
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 50
#> Unobserved stochastic nodes: 2
#> Total graph size: 59
#>
#> Initializing model
jags_fit
#> Inference for Bugs model at "modelo_normal.bugs", fit using jags,
#> 1 chains, each with 10000 iterations (first 1000 discarded)
#> n.sims = 9000 iterations saved
#> mu.vect sd.vect 2.5% 25% 50% 75% 97.5%
#> mu 1.91 0.313 1.31 1.70 1.91 2.12 2.53
#> sigma2 4.75 1.449 3.25 4.07 4.62 5.28 6.91
#> deviance 223.50 2.388 221.46 222.00 222.82 224.27 229.28
#>
#> DIC info (using the rule, pD = var(deviance)/2)
#> pD = 2.9 and DIC = 226.3
#> DIC is an estimate of expected predictive error (lower deviance is better).
# podemos ver las cadenas
traceplot(jags_fit, varname = c("mu", "sigma2"))
Cuando generamos una muestra de la distribución posterior usando MCMC, buscamos que:
Los valores simulados sean representativos de la distribución posterior, esto implica que no deben estar influenciados por el valor inicial (arbitrario) y deben explorar todo el rango de la posterior.
Debemos tener suficientes simulaciones de tal manera que las estimaciones sean precisas y estables.
Queremos tener un método eficiente para generar las simulaciones.
En la práctica, intentamos cumplir lo más posible estos objetivos, pues aunque en principio los métodos MCMC garantizan que cadenas de Markov infinitamente largas lograran una representación perfecta de la distribución posterior (convergen), siempre debemos tener un criterio para cortar la cadena (finita) y evaluar la calidad de las simulaciones (velocidad de convergencia).
Para determinar la convergencia de la cadena de Markov es conveniente realizar más de una cadena; buscamos ver si realmente se ha olvidado el estado inicial y revisar que algunas cadenas no hayan quedado atoradas en regiones inusuales del espacio de parámetros.
Consideremos 3 cadenas de Markov:
jags_fit <- jags(
model.file = "modelo_normal.bugs", # modelo de JAGS
inits = rerun(3, jags_inits()), # valores iniciales
data = list(x = x, N = N), # lista con los datos
parameters.to.save = c("mu", "nu", "sigma2"), # parámetros por guardar
n.chains = 3, # número de cadenas
n.iter = 50, # número de pasos
n.burnin = 0, # calentamiento de la cadena
n.thin = 1
)
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 50
#> Unobserved stochastic nodes: 2
#> Total graph size: 59
#>
#> Initializing model
# podemos ver las cadenas
traceplot(jags_fit, varname = c("mu", "sigma2"))
Las gráficas anteriores nos pueden ayudar, mediante una inspección visual, a determinar si elegimos un periodo de calentamiento adecuado o si alguna(s) cadena(s) está(n) alejada(s) del resto.
Además de realizar gráficas podemos usar la medida de convergencia \(\hat{R}\) que la función jags calcula por omisión para una cadena de Markov:
jags_fit$BUGSoutput$summary
#> mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
#> deviance 223.892 4.6145 221.455 222.120 223.033 224.224 228.773 1.18 130
#> mu 1.972 0.3392 1.351 1.772 2.008 2.176 2.497 1.06 150
#> nu 0.219 0.0456 0.132 0.186 0.218 0.249 0.303 1.05 96
#> sigma2 4.902 2.0805 3.295 4.018 4.586 5.378 7.566 1.05 96
La medida \(\hat{R}\) se conoce como el factor de reducción potencial de escala o diagnóstico de convergencia de Gelman-Rubin, éste es la posible reducción en la longitud de un intervalo de confianza si las simulaciones continuaran infinitamente.
Una regla usual es iterar hasta alcanzar un valor \(\hat{R} \leq 1.1\) para todos los parámetros.
Una vez que tenemos una muestra representativa de la distribución posterior, nuestro objetivo es asegurarnos de que la muestra es lo suficientemente grande para producir estimaciones estables y precisas de la distribución posterior.
Para ello usaremos la salida n.eff de la función jags, que es el tamaño efectivo de muestra. Si las simulaciones fueran independientes, n.eff sería el número total de simulaciones; sin embargo, las simulaciones de MCMC suelen estar correlacionadas, el tamaño efectivo nos dice qué tamaño de muestra de observaciones independientes nos daría la misma información que las simulaciones de la cadena de Markov.
Usualmente nos gustaría obtener un tamaño efectivo de al menos \(100\).
Hay varias maneras para mejorar la eficiencia de un proceso MCMC:
Paralelizar: no disminuimos el número de pasos en las simulaciones pero podemos disminuir el tiempo que tarda en correr.
Cambiar la parametrización del modelo o transformar los datos. En el caso de variables continuas suele ser conveniente centrar los datos. Otra opción es utilizar parametrizaciones redundantes.
Adelgazar la muestra: esto nos puede ayudar a disminuir el uso de memoria, consiste en guardar únicamente los \(k\)-ésimos pasos de la cadena. Esto resulta en cadenas con menos autocorrelación.
Gelman recomienda los siguientes pasos cuando uno esta simulando de la posterior:
Cuando definimos un modelo por primera vez, establecer un valor bajo para n.iter, por ejemplo 10 o 50. La razón es que la mayor parte de las veces los modelos no funcionan a la primera por lo que sería pérdida de tiempo dejarlo correr mucho tiempo antes de descubrir el problema.
Si las simulaciones no han alcanzado convergencia aumentamos las iteraciones a 500 o 1000, de tal forma que las corridas tarden segundos o unos cuantos minutos.
Si JAGS tarda más que unos cuantos minutos (para problemas del tamaño que estamos viendo en la clase) y aún así no alcanza convergencia, entonces juega un poco con el modelo (por ejemplo intenta transformaciones lineales). Gelman sugiere más técnicas para acelerar la convergencia en el capitulo 19 del libro Data Analysis Using Regression and Multilevel/Hierarchical models.
Otra técnica conveniente cuando se trabaja con bases de datos grandes (sobretodo para la parte exploratoria) es trabajar con un subconjunto de los datos, quizás la mitad o una quinta parte.
Understanding Computational Bayesian Statistics, William M. Bolstad.
Doing Bayesian Data Analysis, John K. Kruschke.
Data Analysis Using Regression and Multilevel/Hierarchical models, Andrew Gelman, Jennifer Hill.